#实践课6#####
  
if (!requireNamespace("BiocManager", quietly = TRUE))  
  install.packages("BiocManager")  

# 使用BiocManager安装DESeq2包  
BiocManager::install("DESeq2")  

# 加载DESeq2包  
library(DESeq2)  

# 设置工作目录  
setwd("C:\\Users\\29228\\Desktop\\2024春学期各项汇报\\R语言与生物信息学")  

# 读取mRNA表达数据  
dataDESeq <- read.csv("mRNA_exprSet.csv", row.names = 1, check.names = FALSE) # 第一列是行名  

# 移除在所有样本中表达量均为0的基因  
dataDESeq1 <- dataDESeq[rowSums(dataDESeq) != 0, ]  

# 创建colData，其中包含了样本的分组信息  
# 利用ID从第14和15个字符处提取信息  
  
colData1 <- data.frame(condition = substring(colnames(dataDESeq1), 14, 15))  

# 将"01"和"11"替换为"肿瘤"和"正常"  
colData1$condition <- ifelse(colData1$condition == "01", "肿瘤", ifelse(colData1$condition == "11", "正常", NA))  

# 检查是否所有条件都已经被正确替换  


# 创建DESeqDataSet对象  
dds <- DESeqDataSetFromMatrix(countData = as.matrix(dataDESeq1),  
                              colData = colData1,  
                              design = ~ condition)  

# 运行DESeq  
dds1 <- DESeq(dds)  

# 获取结果  
res <- results(dds1)  

# 写入结果到CSV文件  
write.csv(res, "result.csv", row.names = TRUE)  

# 过滤结果，选择padj < 0.05 且 |log2FoldChange| > 1的基因  
res2 <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)  

# 写入过滤后的结果到CSV文件  
write.csv(res2, "result1.csv", row.names = TRUE)

